JSDM models test

Data description

Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
In this file we have compared the 3 models HMSC,GJAM, JSDM on 21 simulated Dataset.\ We have compared the ability to recover the true interactions used for modelling this datasets, by comparing estimated correlation matrix and matrix representing true interactions.
There is a separate file, with the functions needed to fit each model and study the convergence and more parameters are presented.
## Environment filtering 5 species

JSDM

GJAM

HMSC

Convergence:

Study of interactions

Plot species responses

# to be finished
#fundamental niche
np<-500
nspecies<-5
niche_optima  = seq(2, 98, length.out=nspecies)
niche_breadth = 20

#dataframe for fundamental niches
table_fundamental<-data.frame()

for(i in 1:nspecies){
  tmp<-data.frame(xx=0:100,niche=dnorm(0:100,mean=niche_optima[i],sd=niche_breadth)/dnorm(niche_optima[i],mean=niche_optima[i],sd=niche_breadth),species=i)
  table_fundamental<-rbind(table_fundamental,tmp)
}

xx<-data$env

table_jsdm<-data.frame()
  
for(i in 1:nspecies) {
  tmp<-data.frame(xx=xx,mean=pred_j_mean[,i],q_95=pred_j_95[,i],q_05=pred_j_05[,i],type=rep("jsdm",np),species=rep(i,np))
  table_jsdm<-rbind(table_jsdm,tmp)
}


table_gjam<-data.frame()
  
for(i in 1:nspecies) {
  tmp<-data.frame(xx=xx,mean=pred_gj_mean[,i],q_95=pred_gj_95[,i],q_05=pred_gj_05[,i],type=rep("gjam",np),species=rep(i,np))
  table_gjam<-rbind(table_gjam,tmp)
}


table_hmsc<-data.frame()
  
for(i in 1:nspecies) {
  tmp<-data.frame(xx=xx,mean=pred_hm_mean[,i],q_95=pred_hm_95[,i],q_05=pred_hm_05[,i],type=rep("hmsc",np),species=rep(i,np))
  table_hmsc<-rbind(table_hmsc,tmp)
}
  
#table for predictions
table<-rbind(table_jsdm,table_gjam,table_hmsc)

#table for observations
Y_data = subset(data, select = -env)

table_obs<-data.frame()

for(i in 1:nspecies){
  tmp<-data.frame(xx=xx,obs=Y_data[,i],species=rep(i,np))
  table_obs<-rbind(table_obs,tmp)
}


######## FIRST TYPE OF PLOT
# 1 plot for each species, so one figure with 5 plots, in total 3 figures, one for each model

#jsdm 

for(i in 1:5)
local({
  i<-i
  tmp<-table_jsdm[which(table_jsdm$species==i),]
  tmp_obs<-table_obs[which(table_obs$species==i),]
  tmp_fund<-table_fundamental[which(table_fundamental$species==i),]
  g<<-ggplot()+
    geom_ribbon(aes(x=tmp$xx,ymin=tmp$q_05,ymax=tmp$q_95),alpha=0.5)+
    geom_line(aes(x=tmp$xx,y=tmp$mean,color = "Predicted probability"),lwd=1.5)+
    geom_point(aes(x=tmp_obs$xx,y=tmp_obs$obs),col="#000066",size=0.5) +xlab("Environmental gradient")+ylab("Probability of presence")+
    geom_line(data=tmp_fund,aes(x=tmp_fund$xx,y=tmp_fund$niche,color = "Fundamental niche"),lwd=1)+
    labs(title=paste0("JSDM, Species ",i))+
    scale_color_manual(name = c("Legend"), values = c("Predicted probability" = "#FF6666","Fundamental niche"="#9999FF"))

    assign(paste0("p",i), g, pos =1)
})


grid.arrange(p1,p2,p3,p4,p5,nrow=ceiling(nspecies/2))

#gjam
for(i in 1:5)
local({
  i<-i
  tmp<-table_gjam[which(table_gjam$species==i),]
  tmp_obs<-table_obs[which(table_obs$species==i),]
  tmp_fund<-table_fundamental[which(table_fundamental$species==i),]
  g<<-ggplot()+
    geom_ribbon(aes(x=tmp$xx,ymin=tmp$q_05,ymax=tmp$q_95),alpha=0.5)+
    geom_line(aes(x=tmp$xx,y=tmp$mean,color = "Predicted probability"),lwd=1.5)+
    geom_point(aes(x=tmp_obs$xx,y=tmp_obs$obs),col="#000066",size=0.5) +xlab("Environmental gradient")+ylab("Probability of presence")+
    geom_line(data=tmp_fund,aes(x=tmp_fund$xx,y=tmp_fund$niche,color = "Fundamental niche"),lwd=1)+
    labs(title=paste0("GJAM, Species ",i))+
    scale_color_manual(name = c("Legend"), values = c("Predicted probability" = "#FF6666","Fundamental niche"="#9999FF"))

  
  assign(paste0("p",i), g, pos =1)
})


grid.arrange(p1,p2,p3,p4,p5,nrow=ceiling(nspecies/2))

#hmsc
for(i in 1:5)
local({
  i<-i
  tmp<-table[which(table_hmsc$species==i),]
  tmp_obs<-table_obs[which(table_obs$species==i),]
  tmp_fund<-table_fundamental[which(table_fundamental$species==i),]
  g<<-ggplot()+
    geom_ribbon(aes(x=tmp$xx,ymin=tmp$q_05,ymax=tmp$q_95),alpha=0.5)+
    geom_line(aes(x=tmp$xx,y=tmp$mean,color = "Predicted probability"),lwd=1.5)+
    geom_point(aes(x=tmp_obs$xx,y=tmp_obs$obs),col="#000066",size=0.5) +xlab("Environmental gradient")+ylab("Probability of presence")+
    geom_line(data=tmp_fund,aes(x=tmp_fund$xx,y=tmp_fund$niche,color = "Fundamental niche"),lwd=1)+
    labs(title=paste0("HMSC, Species ",i))+
    scale_color_manual(name = c("Legend"), values = c("Predicted probability" = "#FF6666","Fundamental niche"="#9999FF"))

  assign(paste0("p",i), g, pos =1)
})


grid.arrange(p1,p2,p3,p4,p5,nrow=ceiling(nspecies/2))

########SECOND TYPE OF PLOT
for(i in 1:5)
local({
  i<-i
  tmp<-table[which(table$species==i),]
  tmp_obs<-table_obs[which(table_obs$species==i),]
  tmp_fund<-table_fundamental[which(table_fundamental$species==i),]
  g<<-ggplot()+
    geom_ribbon(aes(x=tmp$xx,ymin=tmp$q_05,ymax=tmp$q_95,col=as.factor(tmp$type)),alpha=0.5)+
    geom_line(aes(x=tmp$xx,y=tmp$mean,col=as.factor(tmp$type)),lwd=1.5)+
    geom_point(aes(x=tmp_obs$xx,y=tmp_obs$obs),col="#000066",size=0.5) +xlab("Environmental gradient")+ylab("Probability of presence")+
    geom_line(data=tmp_fund,aes(x=tmp_fund$xx,y=tmp_fund$niche,color = "Fundamental niche"),lwd=1)+
    labs(title=paste0("HMSC, Species ",i))+
    scale_color_manual(name = c("Legend"), values = c("jsdm" = "#FF6666","gjam" = "#FFFF66","hmsc" = "#FFB266","Fundamental niche"="#9999FF"))
  assign(paste0("p",i), g, pos =1)
})


grid.arrange(p1,p2,p3,p4,p5,nrow=ceiling(nspecies/2))

ALL3<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
 par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
  corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R  JSDM ")
  corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R GJAM")
  corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R HMSC")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(5))

ALL4<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
 par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
  #corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  #title("Partial correlation  JSDM ")
  corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation GJAM")
  corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation HMSC")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(5))

Environmental filtering 10 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(10))

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(10))

Environmental filtering 20 species

JSDM

Gjam

Gjam dimension reduction

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(20))

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(20))

Environmental filtering + Facilitation dense 5 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[4]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(10))

Environmental filtering + Facilitation dense 10 species

JSDM

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[5]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[5]])

Environmental filtering + Facilitation dense 20 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[6]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[6]])

Environmental filtering + Facilitation sparse 5 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[7]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[7]])

Environmental filtering + Facilitation sparse 10 species

JSDM

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[8]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[8]])

Environmental filtering + Facilitation sparse 20 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[9]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[9]])

Environmental filtering + Competition dense 5 species

JSDM

Gjam

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[10]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[10]])

Environmental filtering + Competition dense 10 species

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[11]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[11]])

Environmental filtering + Competition dense 20 species

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[12]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[12]])

Environmental filtering + Competition sparse 5 species

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[13]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[13]])

Environmental filtering + Competition sparse 10 species

JSDM

Gjam

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[14]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[14]])

Environmental filtering + Competition sparse 20 species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[15]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[15]])

Environmental filtering + Competition +Facilitation 5 dense species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])

Environmental filtering + Competition +Facilitation 10 dense species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])

Environmental filtering + Competition +Facilitation 20 dense species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])

Environmental filtering + Competition +Facilitation 5 sparse species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])

Environmental filtering + Competition +Facilitation 10 sparse species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])

Environmental filtering + Competition +Facilitation 20 sparse species

JSDM

Gjam

HMSC

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])